/*
 * Copyright (c)	2010	Bo Zhou<Bo.Schwarzstein@gmail.com>
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 
		http://www.apache.org/licenses/LICENSE-2.0 
 
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

#include "TaubinMeshSmoothing.h"

TaubinMeshSmoothing::TaubinMeshSmoothing(const MObject& MeshObject) : LaplacianMeshSmoothing(MeshObject), mStep(0)
{
}

TaubinMeshSmoothing::~TaubinMeshSmoothing()
{
}

void TaubinMeshSmoothing::SetMu(const float Mu)
{
	mMu = Mu;
}

const float TaubinMeshSmoothing::GetMu() const
{
	return mMu;
}

bool TaubinMeshSmoothing::Smooth()
{
	const float Lambda = GetLambda();
	const float Mu = GetMu();
	const int Step = mStep;

	const WeightType WT = GetWeightType();

	float DiffusionDirection = 0.0f;
	float DiffusionCoeff = 0.0f;

	if( (Step & 1) == 0 )
	{
		DiffusionDirection = 1.0f;
		DiffusionCoeff = GetLambda();
	}else
	{
		DiffusionDirection = -1.0f;
		DiffusionCoeff = GetMu();
	}

	// Do smoothing.
	// 
	const int VertexCount = mMesh.numVertices();
	MPointArray NewPoints;
	mMesh.getPoints(NewPoints);

	std::valarray<int> VertexIds( VertexCount );
	for(size_t i=0; i<VertexIds.size(); ++i)
	{
		VertexIds[i] = (int)i;
	}

	if( IsUsingPaintedVertices() )
	{
		VertexIds = mPaintedVertices;
	}

	for(size_t i=0; i<VertexIds.size(); ++i)
	{
		const int vid = VertexIds[i];
		std::set<int> Neighbors;

		// Get neighborhood of indices
		// 
		int* p = &mVertexSharedPolygons[vid*(mMaxVertexSharedTime+1)];
		while( *p > -1 )
		{
			const int Polygon = *p;
			MIntArray IndicesOnPolygon;
			mMesh.getPolygonVertices(Polygon, IndicesOnPolygon);
			for(unsigned int j=0; j<IndicesOnPolygon.length(); ++j)
			{
				Neighbors.insert( IndicesOnPolygon[j] );
			}
			IndicesOnPolygon.clear();
			IndicesOnPolygon.setLength(0);
			++p;
		}

		std::set<int>::iterator itr = Neighbors.find(vid);
		assert( itr != Neighbors.end() );
		Neighbors.erase(itr);

		std::valarray<double> Weights(1.0/(double)Neighbors.size(), Neighbors.size());

		// Calculate new position.
		//
		MPoint CentralPoint;
		mMesh.getPoint(vid,CentralPoint);
		if( WT == LaplacianMeshSmoothing::EQUAL )
		{
			MPoint WeightedPoint;
			size_t WeightsIndex = 0;

			for(itr=Neighbors.begin(); itr!=Neighbors.end(); ++itr)
			{
				MPoint NeighborPoint;
				mMesh.getPoint(*itr, NeighborPoint);

				WeightedPoint += (NeighborPoint - CentralPoint)*Weights[WeightsIndex];

				++WeightsIndex;
			}

			NewPoints[vid] = CentralPoint + DiffusionDirection*DiffusionCoeff*WeightedPoint;
		}else if( WT == LaplacianMeshSmoothing::FUJIWARA )
		{
			MPoint WeightedPoint;
			size_t WeightsIndex = 0;

			for(itr=Neighbors.begin(); itr!=Neighbors.end(); ++itr)
			{
				MPoint NeighborPoint;
				mMesh.getPoint(*itr, NeighborPoint);

				MVector E = NeighborPoint - CentralPoint;

				double Weight = 1.0 / E.length();
				WeightedPoint += E*Weight;

				Weights[WeightsIndex] = Weight;
				++WeightsIndex;
			}

			WeightedPoint.x /= Weights.sum();
			WeightedPoint.y /= Weights.sum();
			WeightedPoint.z /= Weights.sum();

			NewPoints[vid] = DiffusionDirection*DiffusionCoeff*WeightedPoint;
		}else
		{
			assert(0);
			NewPoints[vid] = CentralPoint;
		}
	}

	mMesh.setPoints(NewPoints);
	
	NewPoints.clear();
	NewPoints.setLength(0);

	++mStep;

	return true;
}
